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1 INTRODUCTION 



Advances in submillimeter telescopy have enabled the 
discovery of flattened structures, in some cases clearly 
Keple rian disks, sur r ounding massive ( ^ IOMq) proton - 
stars dCesaronil 120051 : IChini et al l 12004 IPatel et all 120051 ; 
iBeuther et al. 20061 ). The logical inference - that high mass 
star formation (HMSF) proceeds through disk accretion - 
raises a question: can such disks process rapid mass accre- 
tion, or do they fragment to produce secondary stars? We 
seek to answer this question by estimating the criterion for 
disk fragmentation in the vicinity of a massive star. 

If massive-star disks typically do not fragment, then 
disk accretion poses little barrier to massive star formation. 
Conversely if disk fragmentation occurs, then accretion onto 
th e central star may be (p artially) choked off, as suggested 
bv lTan fc Blackmanl l|2005l ) in the context of low-luminosity 
active galactic nuclei. Moreover each massive star that forms 
from a fragmenting disk may be surrounded by smaller stars 
that began as disk fragments. It is therefore important to 
evaluate models for HMSF in light of the disk fragmentation 
criterion. 

The current work uses and extends the results of 
iMatzner fc Levinl l|2005l . hereafter ML05) , who showed that 
protostellar disks around low mass stars are strongly sta- 
bilized against fragmentation by a combination of viscous 
heating and irradiation by the central protostar. Massive 
star formation is fundamentally different, however, in three 
important ways. The rapid rise of luminosity with mass im- 
plies that stellar irradiation is far more intense; this tends to 
stabilize disks against fragmentation. H owever, massive stars 
must accrete quickly to form at all (e.g.. lWolfire fc CassinellH 



ABSTRACT 

We examine whether massive-star accretion disks are likely to fragment due to self- 
gravity. Rapid accretion and high angular momentum push these disks toward frag- 
mentation, whereas viscous heating and the high protostellar luminosity stabilize 
them. We find that for a broad range of protostar masses and for reasonable accretion 
times, massive disks larger than ~ 150 AU are prone to fragmentation. We develop 
an analytical estimate for the angular momentum of accreted material, extending the 
analysis of Matzner and Levin to account for strongly turbulent initial conditions. In 
a core-collapse model, we predict that disks are marginally prone to fragmentation 
around stars of about four to 15 Mq - even if we adopt conservative estimates of 
the disks' radii and tendency to fragment. More massive stars are progressively more 
likely to fragment, and there is a sharp drop in the stability of disk accretion at the 
very high accretion rates expected above 110 solar masses. Fragmentation may starve 
accretion in massive stars, especially above this limit, and is likely to create swarms 
of small, coplanar companions. 

Il987l) . and rapid accretion favors fragmentation. These ef- 
fects compete to set the critical radius outside of which 
disks fragment. Whether fragmentation actually occurs de- 
pends on the initial disk radius, which itself depends on 
the physical state of the gas prior to accretion. We discuss 
disk fragmentation in 321 considering the stabilizing effect 
of viscous heating ( <I2.2|I before incorporating irradiation by 
the central star ( H2.3p . This combination allows us to iden- 
tify (H2.4C the disk radiu s at which fragmentation sets in. 
The iMcKee fc Tanl (|2003h core collapse model is examined 
in more detail in SJ31 we compute angular momentum scales 
in t|3.2l using formulae derived in the Appendix. In we 
calculate expected fragmentation radii for a range of masses 
in the core collapse model. 

Turning to the consequences of fragmentation, we exam- 
ine in ^5]the likely properties of stars born within fragments 
and the possibility that fragmentation limits accretion. In 
H5.2l we compare our results with observed regions of HMSF. 



2 DISK FRAGMENTATION 
2.1 Criterion for Fragmentation 

We shall concentrate on fragmentation due to local gravita- 
tional instabilities, which set in when Toomre's parameter 

~ c ad f2 
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descends toward unity. Here E is the disk's surface density, 
Q its orbital frequency, and c at j is its adiabatic sound speed. 
(We shall frequently refer to the isothermal sound speed 



7 



-1/2 



Cad where 7 is the ratio of specific heats.) 
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Observatio nal inferences of massive circumstellar tori 
|Cesaronill20"05h indicate that they may be subject to global 
instabilities as well. We discuss this possibility briefly in i )4.4l 
however the local instabilities tend to occur first, and their 
relation to fragmentation is better understood. 

Several authors have identified the frag mentation 
bound ary in terms of a cooling time. Following iGammiel 
l|200ll ) and ML05, we convert this criterion to a critical 
mass accretion rate at a given midplane temperature. The 
cooling time r c is the ratio of the internal energy per area, 
U = c^ d E/[7(7 — 1)], to the dissipation rate per unit area - 
which, in steady accretion, is 



2F V = 



3Affl 2 
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(2) 



where F v is the flux through each disk face. Eliminating E, 
the maximum accretion rate is 



Mn 
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For later convenience we write this in terms of the is othermal 
sound speed c 3 = 7 _1/2 c ad . Extrapolating from iGammid 
<|200ll )'s two-dimensional simulations for a stiff equation of 
state, ML05 estimate M max = 0.89cl/G for the case of a 
three-dimension al, 7 = 5/ 3 disk . Using smoothed-particle 
hydrodynamics, I Rice et all (|2005T l have simulated just such 
a disk, finding Q.t c to lie between 6 and 7 when it frag- 
ments. Assuming also Q = 1, this implies M max ~ (0.37 to 
0.43)Cs/G. For a fixed M, the critical temperature is then 
significantly (f .6 times) higher than ML05 estimated. 

At face value this weakens the conclusion, reached by 
ML05, that fragmentation is unattainable in low-mass pro- 
tostellar disks. However, more recent simulations show that 
disks remain stable at shorter r c depending on how abruptly 
cooling is implemented (E. Harper-Clark 2006, private com- 
munication). Due to uncertainties in the aforementioned 
cooling factor, and in light of the stringent reso lution re- 
quirements for collapse outlined by Nelson! J2006h . w e con- 
sider the value of D.t c obtained bv lRice et al.l |2005) to be 
uncertain by up to a factor of two. We shall therefore be 
conservative, by adopting the stricter ML05 criterion that 
fragmentation occurs when 

Cs < C s , crit = 1 m{GKL d f /3 . (4) 

Since the mass accretion rate is comparable to ec c ff (core) 3 /G 
in the core collapse scenario, where c c ff (core) 2 is the ratio of 
pressure to density in the core, equation @ implies, quali- 
tatively, that a disk fragments if its sound speeds falls below 
the effective sound speed of its parent core (see equation 
|23| below). This point was made for thermally supported 
cores by ML05, and equation simply extends the rule to 
turbulent cores. 

Equation Q is conservative in the sense that fragmen- 
tation may also occur at somewhat higher values of c 3 . It 
is even more conservative given that 7 declines from 5/3 at 
the higher temperatures relevant to massive star formation. 
We shall make several other conservative estimates in or- 
der to show that disk instability is all but inevitable during 
massive star formation. 

When assessing disk stability in a given scenario, 
we first calculate the midplane temperature profile c a (r) 
of the disk given its central mass Af*, central luminos- 
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Stellar irradiation 
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Viscous flux 



ity L+, and accreti on rate M*d- For this we adopt the 
IShakura fc Sunvaevl a parametrization of viscosity, in which 
the steady mass accretion rate is 



M vis 



n 



(5) 



For the choice of critical temperature made in equation ((4} 
this corresponds to a = 0.30Q at the onset of fragmentation. 
We keep a fixed at 0.30 throughout our analysis, as this 
correctly reproduces the fragmentation boundary, although 
this is an overestimate for stable disks. 

Having identified the fragmentation radius as the loca- 
tion where the disk sound speed falls to the critical value 
in equation Q, we then compare this to the characteristic 
disk radius 



R d 



GAL 



(6) 



for accreting gas with specific angular momentum j. 

We set the viscous accretion rate equal to the accre- 
tion rate from the envelope, and let Q — 1 . We then check 
whether the disk can remain in this steady state with two 
models for heat generation. In i|2.2l we ignore the luminosity 
from the protostar and find a minimum value for T<j(r) using 
only the viscous generation of heat. Next, in £|2.3I we include 
the flux from the protostar received at the disk surface, and 
again solve for the midplane temperature as a function of 
radius. 



2.2 Viscous Heating 

In a thermal steady state, the flux of viscous energy radiated 
by each face of the disk is given by equation (01. All of the 
disks considered in this paper are optically thick to their own 
thermal radiation; therefore, the flux can also be derived 
from radiation transfer across an optical depth kE/2 from 
the disk midplane to its surface: 

8 

3t r ~ 



F, 



-crT d , 



(7) 



where tr = KijEd/2 is the optical depth corresponding to 
the Rosseland mean opacity Kn(Td). The factor 8/3 in equa- 
tion (0 is derived by assu ming that the dissipa tion rate 
per unit mass is a constant l|Chick fc Cassen)ll997l). We ob- 
tain t emperature dependent opacities from ISemenov et all 
(2003). These opacities are very insensitive to density; we 
adopt values for 10 -12 ' 5 gem -3 , an appropriate value for a 
disk with Q = 1 and a period of a few centuries. 

Neglecting irradiation of the disk surface, in a steady 
state F v — F r and M v isc = M+d- We solve for E in equation 
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((5|, and then c s from equations ([2]) and ([7]) using pc 2 a = UbT 
(for molecular weight fi), yielding 



3k%K R MLQ 3 



(8) 



1287T 2 <j/i 4 a 

This is an implicit formula for c s , as kr depends on temper- 
ature. Equating c s to c s . cr it (equation [3| implies fragmenta- 
tion when 



fi < 8.54 



G 10/3^4/3^1/3 



(9) 



In practice we calculate the run of c s (r) and kr(t) self- 
consistently in order to evaluate equation JSJ. 

Whereas ML05 found a unique value of f2 C rit for opti- 
cally thick accretion disks around low-mass protostars, we 
shall show below that equation (|9} gives a roughly constant 
fragmentation radius of about 130 AU. The key difference is 
the higher accretion rate, which implies much higher critical 
temperature (hundreds of K) in massive star formation com- 
pared to ~ 16 K for the low mass case. As the opacity law 
does not obey kr(T) oc T 2 for these higher temperatures, 
the ML05 result does not hold for these more massive stars. 



2.3 Stellar Irradiation 

Stars of mass somewhat greater than IOMq undergo Kelvin- 
Helmholtz contraction rapidly enough to settle onto the hy- 
drogen burning main sequence while still accreting. As the 
main sequence luminosity increases rapidly with M*, the 
consequences of stellar luminosity become acute for mas- 
sive stars. One such consequence is the heating of the disk 
midplane due to reprocessed stellar radiation. In principle, 
this stabilizes disks to longer periods (larger radii) than we 
found in equation <(9j by considering only viscous heating. 
We parametrize irradiation through the reprocessing factor 
/ defined as the ratio between the incident flux Fi rl normal 
to the disk surface, and the spherical stellar flux at that 
radius: 



F — 

i irr — 



(10) 



A calculation of / requires a model for the reprocessing of 
starlight onto the disk. 



2.3.1 Infall geometry 

ML05 model / in the context of an infall ge o metry sim- 
ilar to that used by IWhitnev fc Hartmannl l|l993l ) and 
iKenvon et ai1l|l993h for scattered-light images of protostars. 
Starlight is absorbed at the inner edge of a rotating infall 
envelope l|Chevalierll 19831 ; iTerebev et al.lll984h whose inner- 
most streamlines have been removed by the ram pressure of 
a magnetically collimated protostellar wind. This removal 
implies a suppression of the star-disk accretion rate relative 
to the rate at which mass would otherwise accrete - from the 
surrounding core, for instance, in a core-accretion model: 



Ah 



eM c ; 



(11) 



where e = cos 9o if streamlines are removed from all angles 
within 6*o of the axis. (We assume the prestellar matter is 
isotropically distributed before it falls in.) 



The balance of forces that de termines e is modeled in 
detail bv lMatzner fc McKed {2000) for low-mass star forma- 
tion, and we expect that their analysis holds into the massive 
star regime. The location of the innermost streamline is a 
simple function of e: it strikes the disk at a radius 



(l-e 2 )^. 



(12) 



The calculation of / is particularly simple if the dust 
envelope is (1) optically thick to stellar photons; (2) optically 
thin to dust thermal radiation; and (3) not hot enough to 
sublimate. Under these conditions, ML05 find 



/ ~ O.le" ' 35 

for a reasonable range of e. 



(13) 



2.3.2 Envelope self-opacity and dust sublimation 

It is important to examine the assumptions that led to equa- 
tion (|13l) . especially in the context of very rapid accretion 
onto very massive stars. 

In ij3.2l we will estimate the typical infall column den- 
sity S sp h,/ in a core collapse scenario, and find it to be a 
few times lower than t he core column E e , i.e ., £ sp h,/ — 0.3 
g cm - . Given that the Scmcnov et all (|2003l ) Planck opac- 
ity peaks between 3 and 16 cm 2 g _1 for 45 < T < 930 K, we 
expect the disk midplane to be shielded by moderate optical 
depths (~ 1 — 5) from the reprocessing surface. A solution to 
the radiation diffusion problem is beyond the scope of this 
paper. Instead we will apply equation (|13[) ; this is conserva- 
tive, in the sense that it overestimates the stabilizing effect 
of irradiation. 

A further potential complication arises for especially 
high accretion rates, those in excess of 1.7 x lO~ 3 Af0 yr _1 . 
The critical temperature of such a disk is quite hot, T cr i t > 
1050 K according to equation Q. As we shall see in f|4l 
the stabilizing effect of irradiation is accentuated in these 
disks by a sharp drop in Rosseland opacity. If the disk 
is to reach 1050 K, however, dust in the nearby infall 
envelope must approach the silicate sublimation tempera- 
ture of about 1500 K. This leads to the disappearance of 
dust within th e infall envelope inside a certa in sublima- 
tion radius R 3 . iMonnier fc Millan-Gabetl <|2002ft determine 
Rs ~ 35(L*/1O 6 L ) 1/2 AU by optical interferometry (and 
note that this value requires the existence of large silicate 
grains). We expect that disks with Rd > Rs are relatively 
unaffected by dust sublimation. We do not attempt to cal- 
culate disk irradiation in cases where this is not true. 



2.3.3 Incorporation into fragmentation calculation 

Solving equation (|10p . we can find the equilibrium temper- 
ature of the outer reaches of an optically thick disk for 
which irradiation dominates over viscous heating. While 
low-mass protostellar disks can be stabilized in this regime 
(ML05), massive-star disks are not. We therefore account 
self-consistently for Td(Rd), in two steps. The disk's effec- 
tive (surface) temperature (T 3 , say) is determined by the 
requirement that it emit the viscous flux in addition to re- 
emitting the incident flux: 



oTt = F v + F ir 



(14) 
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The midplane temperature Td is derived from T s from radi- 
ation diffusion of the viscous flux across optical depth tr, 

<?T d {R d ) 4 = aT* + ^t r F v = (jtr + l\F v + F irr 

^ ^T R F V +F il r. (15) 

We account for the temperature dependence of the opacity 
when solving this equation numerically. 

As before, we identify the critical disk radius i? cr it at 
which Td(Rd) = Tcrit- Fragmentation occurs if the disk ex- 
tends beyond -Rcrit. 

2.3.4 Other Considerations 

We pause to address two minor concerns: 
Shock Heating- Infalling matter is decelerated in an accre- 
tion shock upon reaching the disk, and heat radiated by this 
shock warms the disk surface. However the gravitational po- 
tential at Rd is very small compared to that at the stellar 
surface, as Rd S> R*. Moreover, the star's emission is domi- 
nated by hydrogen burning rather than accretion, and a fair 
fraction of the starlight is reradiated onto the disk surface in 
our model (eg. I13p . Shock heating is thus wholly negligible 
(by about four orders of magnitude, for a 30 Mq star). 
Radiation Pressure- When Q = 1 the ratio of gas to radia- 
tion pressure at a characteristic fragmentation radius is 

Pg _ 3k B Q 2 = 1.8 30 Mq / tacc V / 150AU \ 3 
P rad ~~ 27vGa^ Tit ~ M* ^ 10 5 yr J \ R d J ' 

where i acc is the duration of accretion (see § H2.4I and 13. 
Radiation pressure remains negligible out to periods of 3300 
years (scaling as (M*E)~ 3//4 in the core model of Sj3|. More- 
over the photon diffusion time across the scale height H, 
fdift — 3trH/c, is a few hundred times shorter than the or- 
bital period. Consequently photon pressure is irrelevant for 
disk fragmentation during massive star formation. 

2.4 Typical Parameters for Fragmentation 

Before treating the core accretion model in detail in fJ21 we 
wish to draw a few conclusions that are reasonably indepen- 
dent of a scenario for massive star formation. We adopt in 
our irradiation model a fiducial efficiency parameter e = 0.5, 
but we consider other values in i|4.2l 

We begin by mapping the maximum disk radius and 
maximum disk angular momentum as functions of the stel- 
lar mass and the accretion time, t acc = 2M*/M t< !. (The 
factor of two derive s from the core accretion scenario of 
iMcKee fc Tanll2003l ; accretion time is simply a convenient 
parametrization for accretion rate.) The results for R CY it and 
jcrit are shown as the solid curves in figures [T] and [2] respec- 
tively. 

Not all of this parameter space is relevant, however. 
Observations of protostellar outflows emerging from sites of 
massive star formation imply dynamical ages of order 10 
years. On both of these figures, we highlight within the dot- 
ted lines a plausible range of t a cc as the range of values 
predicted bv lMcKee fe~T an (2003) for core column densities 
(E c ) of 0.3 - 3 g cm -2 . One may further restrict one's at- 
tention to masses between 10 and 120 Mq, as more massive 




10 100 1000 

M, (Af Q ) 



Figure 1. The fragmentation radius (in AU) is shown for a range 
of central star masses and accretion times. Also shown is the star 
formation time in the core model [dotted lines, for 0.3 < S c i < 3 
g cm -2 ) and the region affected by dust sublimation in the infall 
envelope (filled region), where our model does not hold. The sharp 
kink in the lines of constant radius is due to a drop in opacity for 
disk temperatures above about 1050K. 



# in ng. m 


i°gioi ( cm2 s x ) 


Reference 


1 


22.3 


Cesaroni et al. (20051 


2 


21.4 


Patel et al. (2005) 


5 


22.1 


Zhang et al. (2002) 


6 


22.0 


Bernard et al. ( 1999) 


9 


21.6 


De Buizer & Minier (2005) 



Table 1. Estimates of angular momentum of observed disks in 
figure [3] and corresponding references. Angular momentum esti- 
mates are computed from the observed velocity gradient over the 
extent of the disk. Estimates are only made for data points that 
showed a clear velocity gradient associated with a disk or torus. 



stars are not known to exist. From this we infer: for plausi- 
ble values of the accretion time scale, massive-star accretion 
disks fragment for radii above about 100-200 AU, with 150 
AU being a typical threshold. 

Table ([TJ lists our estimates for the angular momentum 
of several observed disks, most of which appear likely to 
fragment. 

3 CORE ACCRETION 

With these relatively model-independent results in hand, 
let us evaluate disk fragmentation in the turbul ent core 
model (|Mvers fc Fullerl 1 19921 : IMcKee fc Tanl l2003h . which 
posits that massive stars accrete from hydrostatic structures 
(cores, subscript "c") that have assembled from an over- 
dense region (clump, subscript "cl") within a larger molecu- 
lar cloud. This is similar to the standard model of low-mass 
star formation, except that turbulent motions are subsonic 
within low- mass cores (e.g., ML05) and supersonic in mas- 
sive cores. 

To calculate core and core collapse properties, we adopt 
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Conservative^ 
**55^ore, model 
^gtable^ 



22.25. 




a{r) 



Pd = 0.88GE 2 



clumps that host massive star f ormation llPlume et al. 19971; 



Hille nbrand fc Hartmann 1998; Figer et al.lll999l;lKim et all 



200d; Ivan den Bergh et"aH ll99lT[Gilbert fc Graham! l200ll: 



Rc = 



0.040 



•V2£ 



1/2 



m. 



30M F 



1/2 



pc 



1000 



Figure 2. Angular momentum and disk fragmentation. The crit- 
ical value of angular momentum, labeled as log 1 Q(j cr it), is plotted 
as bold curves for a range of stellar masses and formation times, 
tacc = 2M* / (M +d)- Dashed lines show our predicted angular mo- 
mentum for the iMcKee fc Tanl (2003) core model (e = 0.5, k p = 
1.5). As discussed in £13.21 this is also an approximate upper limit 
to j given and t aC c ■ Disks tend to fragment except in the up- 
per, dark filled region. In the lower, light filled region, envelope 
dust grains sublimate within R^; our model is not secure here. 
Dotted lines delimit core model formation times for 0.3 < E c ] < 3 
g cm -2 . 



the IMcKee fc Tanl (|200St ) models and fiducial parameters. 
Turbulent cores are modeled as singular polytropic spheres 
with density profiles p(r) oc r~ kp , with k p ~ 1.5. Since the y 
are assumed to be hydrostatic (see also iTan et al.l 12006), 
they must be pressu re confined within their parent clumps. 
IMcKee fc Tarj|2003l evaluate the mean hydrostatic pressure 
within clumps to be 



(17) 



Further, they assume that massive cores are segregated to 
the centers of clumps where the pressure is roughly twice 
Pel- Fixing cores' surface pressure to this value implies, 
for fiducial parameters, that E c = 1.22E c i. Observations 
imply a column density E c j of around 1 g cm" 2 for the 



de Marchi et alj 1 19971 ; iTurner et all |2000| ; iFaundez et all 
20041 ). 

As cores a re bound and argued not to fragment 
l|Krumhol j I2006T ) . their gas either accretes on the star (a 
fraction e) or is blown out by the protostellar wind (1 — s). 
The final stellar mass is therefore M+f = sM c , with e ~ 0.5, 
and the core radius is 



(18) 



where E c i, C gs is in g cm - . 

In the following sections we shall refer to a = cr(r), the 
local velocity dispersion of core gas at radius r; in a singular 
polytropic core, 



47T 

&4>B{k p - 



GM c (r) 



1) 



(19) 



where M c (r) is the enclosed mass, and 4>b = P/(po~ ) 
accounts for magnet ic contributions to the total pressure; 
IMcKee fc Tanl (120031 ) estimate (f> B ^ 2.8. 

We note, in passing, that the fiducial core model im- 
plicitly assumes that stellar accretion halts due to a lim- 
ited mass supply, rather than due to the onset of vigorous 
stellar feedback, fragmentation of the core or disk, or any 
other dynamical effect. Alternatively one could either have 
e = M+f /M c <C 1, or one could define the core to be the sub- 
region that successfully accretes (for which e ~ 0.5). In the 
latter choice, the pressure-equilibrium column, ~ 1.22E c x, 
would provide only a lower limit to E c . As it is widely held 
that the upper mass cutoff derives from stellar feedback, 
we expect columns in excess of this lower limit to prevail 
for very massive stars. Nevertheless, we evaluate the above 
equations for E c = 1.22 gem" 2 in the fiducial case. 



3.1 Collapse 

IMcKee fc Tanl i|2003l ) show that the accretion time is close 
to the free-fall time evaluated at the core's surface density. 
Their equations (3), (4), (5), (35), and (36) specify 



M*d = £</>acc 77 
Gr 



(20) 



where 



0.71(3.38 - k p ){k p 

x 



(l+Ho)V2> 

1, 



-l) 3/2 (3-M 

magnetic 
nonmag. 



1/2 



1.9 



(21) 



where the magnetic term includes levitation by the static 
field (the factor (1 + Ho) — 2) as well as turbulent magnetic 
pressure (the factor 4>b)- The second line evaluates the first 
for a magnetized core with k p = 1.5, for which the final 
accretion time is 



-3/4 



cl,c 



0.5 



1/4 



e 3OM / yr ^ 

jMcKee fc Tanll2003l ). Equation (|21[) implies fragmentation 
when 



< l.O4(e0 



, 1/3 a 

acc / u •> 



(23) 



i.e., for c s < 1.02cr when e = 0.5 and acc = 1.9. 

During collapse, distribution of mass, M c (r) oc r 3 ~ kp 
and free-fall time, ts(r) oc p(r) _1/ ^ 2 , imply 



(t/t*/)*" M* f 
oc M 



and 



t 



where 



(24) 



The accretion time is longer than the Kelvin-Helmholtz 
time for stars > 10M Q (|Wolfire fc Cassinellil [19871 ). imply- 
ing that massive protostars reach the zero age main-sequence 
(ZAMS) during accretion. When calculating the time evolu- 
tion of fragmentation in £14.31 we employ the lMcKee fc Tanl 
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(2003) models for the luminosity of a massive accreting pro- 
tostar; to calculate irradiation at the end of accre tion we also 
includ e ZAMS luminosity, using formulae from iTout et al.l 
<| 19961 ). 



3.2 Angular Momentum 

To evaluate our fragmentation criterion, we require the an- 
gular momentum of material within a colla psing turbulent 
core. Although this is not included in the iMcKee fc Tanl 
l|2003h models, we shall estimate angular momentum within 
these models by generalizing a calculation by ML05 (itself 
an analytical version of the iBurkert fc Bodenheimerl |2000| 
calculation, in the vein of iFleckll 19871 ) . This calculation (pre- 
sented in the Appendix) notes that a core model specifies 
the turbulent velocity dispersion a(r) as well as the density 
profile p(r). For a special class of velocity fields one can com- 
pute the ensemble-averaged specific angular momentum and 
velocity dispersion. The velocity field must be isotropic; it 
must have a velocity difference between any two points that 
scales as a ( | r i — r2 1 ) oc | r i — r-2 1 ^ regardless of the underlying 
density distribution; and its Cartesian components must be 
uncorrelated (i.e., transport no average shear stress). The 
ensemble-averaged specific angular momentum and velocity 
dispersion are calculated in equations (I44ll - (|47| l. We define 
the spin parameter 8j for a turbulent region of size R: 



if) 



1/2 



R°{R) R (a(R) 2 } 



1/2- 



(25) 



The rightmost expression involves root-mean-square ensem- 
ble averages over the turbulent velocity field, which are cal- 
culated in the Appendix. The factor fj accounts for the dif- 
ference between the ratio of rms averages and the ratio of 
amplitudes; since the amplitudes are random variables, this 
includes both an overall offset and a disperson. ML05 es- 
timate log 10 fj — — 0.088lg'49 on the basis of a Gaussian 
model for the velocity field. 

Our evaluation of the spin parameter is presented in 
the Appendix and summarized in figure [7] and table [7] For 
a turbulence supported region with p oc r~ kp one must have 



a(r) 2 oc GM(r)/r oc r l 



so 



- k p /2 -> 1/4 (for 

k p -> 3/2). In the IMcKee fc Tanl core collapse model, each 
shell of matter accretes in sequence. For each shell, the spin 
parameter is 



,(shell) 



0.85/3°' 42 /j 



0.47/, 



(26) 



(see equation I49f) . However, the last shell to accrete natu- 
rally has the highest angular momentum; as the disk accu- 
mulates vector angular momentum [26] may overestimate the 
disk radius. Moreover, although collapse is generally super- 
Alfvenic, magnetic braking may sap j somewhat. A rough 
lower bound on 8j is given by the specific angular momen- 
tum of the entire core, which corresponds to 



(27) 



we must decide 



0j (entire core) = 0.50/3°' 55 /j 0.23/j 

(also equation 25J . 

Given upper and lower limits for 
which value to adopt. The remainder of this paper is in- 
tended to establish that fragmentation is inevitable in the 
outer reaches of massive-star disks, so we shall adopt the 
more conservative estimate (eq. I27|l . Please bear in mind 



that the upper bound to Rd is about four times larger. The 
protostellar outflow will tend to remove \aw-j material (see 
also ML05) , but we expect this effect to be rather minor and 
do not evaluate it. 

For convenience we shall define a second rotation- 
related quantity, cj>j, by 



GhL 



(28) 



The factor (j>j defined here is related to the spin parameter 
0j by fa = Ra(R) 2 9j /[GM(R)}. Hydrostatic equilibrium re- 
quires Ra(R) 2 /[GM(R)} = l/[2(k p - i)<j> B ]: therefore 



2{k p - 1)</>E 



0.360, 



(29) 



giving 



0.067 in the fiducial case. 



The disk acquires a final radius which is about l/40th 
of the core radius: 



Rd 



J = 



300 



-R c 

e 



3/2 



f 



30 M Ed,, 



1/2 



AU 



(30) 



in the fiducial, conservative case given by equation (|27p . 
During accretion, the disk radius remains proportional to 
the current radius of accretion; therefore 



Rd 



R, 



--►2/3 



(31) 



at least on average. 

It is useful to know the column density scale in the 
infall for reference in the calculation of disk irradiation. As 
discussed in ML05, this is characterized by £ sp h(.Rd): the 
column outward from Rd in a nonrotating infall of the same 
accretion rate. We find 



-'sph 



(Rd 



2 5 /2( fcp _ 1)0 S 0/ 
0.64E C 



(32) 



i.e., that the final infall column is comparable to the core 
column (see also ML05). Over the course of accretion, 



T, sp h(Rd) 

S ap h(-Rd,/) 



— 1/3 



(33) 



Finally, a note on the relation between a core's den- 
sity profile and its angular momentum scale. Considering 
the range of values 1 < k p < 2, the angular momentum of 
a turbulent sphere decreases sharply toward zero as k p ap- 
proaches 2 and approaches zero, as shown in figure[7] This 
trend is easily understood: when = 0, the velocity differ- 
ence between two points is independent of their separation 
and must therefore be contained in very small-scale motions 
as in an isothermal gas. Indeed, the three-dimensional power 
spectrum scales as k~ 3 ~ 2/3 and contains a divergent energy 
at small scales as — > 0. Angular momentum is dominated 
by the largest-scale motions that fit within the region of in- 
terest, and therefore vanishes if turbulent energy appears 
only at small scales. However, a core whose hydrostatic sup- 
port is effectively isothermal (0 = 0) may still contain an- 
gular momentum due to background turbulence if it is con- 
fined in a turbulent region with > 0. This situation holds 
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for thermally supported cores within turbulent molecular 
clouds, and formed the basis for the ML05 estimate of disk 
radii in low-mass star formation. 

Background turbulence may increase j for turbulent 
cores, as well, if the density profile flattens and the effective 
value of j3 increases across the core boundary. Our calcula- 
tion in the Appendix assumes that velocity field is described 
by the same j3 everywhere, so the corrected value of 8j should 
be intermediate between the core's (3 and that of the parent 
clump. 

3.3 Observations of Rotation in HMSF 



iGoodman et al.l (| 1993f ) observe velocity gradients of dense 
cores, including both low and high-mass cores, using C ls O, 
NH3 and CS as tracers. They find that the ratio of rota- 
tional to gravitational energy, which we call j3 TO t (to avoid 
confusion with a oc r 13 ), takes a rather broad distribution 



around a typical value of 0.02, i.e., log 10 /3 rot 



-1.7. In our 



model for rotation within singular polytropic cores, 



Prot — 



(5 - 2k p f 



8(3-fc p )(fc p -l) 4>i 



(34) 



Using 8j 



0.23/j as we estimated for the entire core 



(eq. [27]) this gives log 10 f3 TO t ^ —2.1 ± 0.7, whereas us- 
ing 6j — 0.47/j as appropriate to the outer shell (eq. [26]), 
log in Piot — — 1.7±0 .7. The agreement would thus be good if 
the lGoodman et al.l observations traced the outermost core 
gas, but this seems unlikely. A better explanation is that 
the observed cores have somewhat flatter density profiles; 
the discrepancy is removed if k p — 1.35 rather than 1.5 (us- 
ing eqs. [42] and [49]). As noted above, embedding the core 
within a clump medium that has a flatter density profile 
(and higher (3) would have the same effect. Protostellar out- 
flows also raise j slightly by removing material on axis. Of 
course, our model for core angular momentum is based on 
an idealized turbulent velocity spectrum, and undoubtedly 
involves some error. 

We include observations of disks and toroids in massive 
star forming regions in figure (J5} and table [5] However it 
should be noted that such comparisons are limited due to 
uncertainty in both the values of quoted parameters and in 
the actual phenomena being observed. Only recently has it 
been possible to achieve the resolution and sensitivity re- 
quired to constrain models of massive star formation. Many 
uncertainties remain when distinguishing between infall, ro- 
tation (Keplerian or otherwise), and outflow. Many of the 
objects show velocity gradients that are consistent with Ke- 
plerian rotation, but could also be attributed to another 
bulk motion, such as infall. Furthermore even when rota- 
tion is indeed Keplerian, it is difficult to distinguish the disk 
edge from its natal co re ijDe Buizer fc M inicr 20051: ICesaronil 
l2005l;|Pe Buizerl|2006ri . This, and confusion with the outflow 
i De Buizerll200a ). may cause disk masses and extents to be 
overestimated. The mass, luminosity and multiplicity of the 
central object involve further uncertainties. 

We have already noted that the angular momentum 
seen in observations coincides with our estimated upper 
bound from the core collapse model. Although we think 
this agreement is likely to be real, it is also possible that 
the observed rotation is about a stellar group rather than 
a single object. This alternative is strengthened by the fact 



that many rotating envelopes are inferred to be more mas- 
sive than any single central object, on the basis of the cen- 
tral luminosity. Such overweight disks are subject to strong 
global self-gravitational instabilities quite distinct from the 
local instability we addressed in Sj2] We refer the reader to 
IShu et all (jl990l ) on this point; see also the discussion in < M.4I 
and that given in ML05. 

It is notable that the three observations most confi- 
dentl y interpreted as thin mas s ive-star disks ( Patel et alJ 



120051 : IDe Buizer fc Minierl 120051 : IShepherd et~"ai1 l200lf ) are 



the best match to our predictions. The other, more extended 
massive disks or toroids are of uncertai n mass and radiu s and 
may enclose many s tars. As noted by ICesaron 3 (|2005l l and 
iBeuther et ail (2006), in these instances we might be seeing 
proto-cluster, rather than protostellar, disks and toroids. 



4 FRAGMENTATION OF CORE-COLLAPSE 
DISKS 

A numerical evaluation of fragmentation is presented below 
in i)4.1l but first we calculate the scalings that govern these 
results. Given the fragmentation criterion from i )2.1l and the 
fiducial core-collapse model described in $3] we can estimate 
the ratio Td/T CI i t by ignoring either irradiation ("active" 
disks, as in 32. 2[) or viscous heating ("passive" disks, using 
323J: 



TdjRd 

TcTlt 



0.15 



0.35 



/ £ 0.5 K H,cg 3 2j cg 3 \ 
V ^.30 J 



1/20 



'0.5 



M3 V 1 '* 
M */,30 2j «g= 



1/4 



(active disk) 
(passive disk) 



(35) 



where = 10 5 L 5 L e , e = 0.5£o.5, «i? = KR, cgs cm 2 g" 1 , and 
E = E C gs gem -2 . 

Note that increasing M*/ destabilizes disks in both 
regimes. High values of E enhance viscous heating relative 
to irradiation, but neither is sufficient to prevent fragmen- 
tation around a 30 Mq star when E ~ lgcm" 2 . However, 
if we take L — 20 and E = 0.034 gem" 2 - values typi- 
cal of the low-mass star formation studied by ML05 - then 
T d (R d ) > T crit for M > 1.3 Mq according to equation (f33|) . 
The evaluation used here neglects several effects treated in 
that paper, such as thermal core support. Nevertheless these 
scalings explain why disks are intrinsically less stable during 
massive star formation than in the low-mass case, as seen in 
detail below. 



4.1 Results 

Fragmentation is expected when Rdj > Rd,ctit> i.e., when 
the disk extends past the critical fragmentation radius at 
some point during formation. Figure [3] compares these two 
radii for a range of stellar masses, using the fiducial, con- 
servative core collapse model (k p — 3/2, e = 0.5, E c i = 
lgcm -2 , 0j = 0.23 fj) to compute Rd- 

To address the lowest-mass stars whose disks can frag- 
ment, we must account for the thermal component of a 
core's hydrostatic pressure and for the accretion luminos- 
ity, both of which are negligible in very massive stars. We 
have (1) included a thermal component (at 20 K) in the 
effective temperature that sets the accretion rate, so that 
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M*d > 10 _5 ' 3 £Mq yr _1 for all masses; (2) added accretion 
luminosity to the ZAMS luminosit y when estimati n g disk 
irradiation; and (3) employed the IPalla fc Stahlerl (|l992T ) 
models (for M+d = 1O~ 4 A40 yr" 1 , which is appropriate) 
to estimate the radius of the accreting protostar. Although 
rather approximate, these amendments are of diminishing 
importance as M*f increases beyond ~ 20Mq. 

Given the expected range of disk radii, all the disks 
presented in figure [3] are candidates for fragmentation. The 
expected disk radius crosses the fragmentation boundary 
for M+f ~ 3.5Mq, and the two remain almost equal un- 
til M*f ~ IOMq; fragmentation is marginal in this range. 
Fragmentation becomes increasingly likely as the mass in- 
creases, though slowly: i? cr it is within a factor of 2 of Rd 
for M*f < 23M Q . For Af*/ > 57M , _R crit drops below the 
range of disk radii implied by the dispersion of fj given 
below equation (|25p - as indicated by the gray region in 
figure [3] The specific masses quoted depend on our model 
for angular momentum, particularly as the critical radius is 
relatively constant in the range 100-150 AU. 

Recall, however, that the disk angular momentum de- 
rives from a turbulent velocity field and is therefore quite 
stochastic. The spread in j predicted by a Gaussian model 
for the velocity field allows for the frequent formation of 
disks twice as large as predicted in equation (|30p . Likewise, 
much smaller disks (by about a factor of nine) can form 
equally easily from a chance cancellation within the core ve- 
locity field. This dispersion in expected radii is indicated as 
a shaded band in figure [3] Remember also that we adopted a 
conservative estimate of the disk angular momentum; other- 
wise, disk fragmentation would have been even more preva- 
lent. Taking these points into account, we can draw a few 
conclusions with relative certainty. 



(i) A significant portion of the O star and early B star pro- 
tostellar disks predicted in the core collapse model are prone 
to fragmentation, although the exact fraction is sensitive to 
the (uncertain) angular momentum scale and fragmentation 
criterion; 

(ii) The tendency of disks to fragment increases with stel- 
lar mass (A'hf); it decreases with higher column densities 
(Ed), and with steeper initial density profiles (k p ). 

(iii) Disks accreting more rapidly than about 1.7 x 
lO _3 M0yr _1 are destabilized by the sharp drop in dust 
opacity at ~ 1050 K, according to equation (j4|). In the fidu- 
cial core model, this occurs above about 110 Mq, close to 
the observed upper limit of stellar masses. More generally, 
this occurs for M*/ > 87(tacc/lO 5 yr)M . 

(iv) At somewhat higher accretion rates, however, dust 
sublimation invalidates our model for starlight reprocessing 
in the infall envelope. 

(v) So long as disks remain optically thick, any effect 
that decreases the Rosseland opacity is destabilizing. For 
instance, low-metallicity disks are less stable than those of 
solar composition. ( Primordial disks are h owever opaque in 
their inner portion s. iTan fe McKe"e]l2004) By adopting the 
(relatively opaque) ISemenov et al.l 1 20031 ) dust opacities, we 
have underestimated disk fragmentation. 
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Figure 3. Relevant radii. The characteristic disk radius predicted 
by (fiducial) core accretion is accompanied by a shaded band il- 
lustrating our (Maxwellian) model for its dispersion. The largest 
stable radius is plotted for comparison; this is accompanied by the 
contributions from pure irradiation (no viscous heat) and pure 
viscosity (no irradiation). The turnover of R cr it above ~ 20Mq 
is due to the scaling of ZAMS mass and luminosity. The refer- 
ences for the observational points are listed in table (|2j. Circles 
represent objects that may best be described as cores, whereas 
diamonds represent those objects whose disks are well resolved. 
Squares indicate objects for which it is unclear whether they are 
rotating, infalling, or both; see i|5.2l for discussion. 



Number in figure [3] 


Reference 


1 


Cesaroni et al. (2005) 


2 


Patel et al. (2005) 


3 


Olmi et al. f2003) 


4 


Olmi et al. (20031 


5 


Zhane et al. (2002) 


6 


Bernard et al. (1999) 


7 


SheDherd et al. (2001") 


8 


Shepherd et al. (2001) 


9 


De Buizer & Minier (2005) 



Table 2. Observational data points in figurc[3]and corresponding 
references. 



4.2 Effect of Varying Efficiency 

Up to this point we have ad opted e = 0.5 as the fiducial 
accretion efficiency, following iMcKee fc Tanl (120031 ). In the 
theory of iMatzner fc McKee l|2000l ). e is set by the ejec- 
tion of material by a centrally-collimated protostellar wind. 
IMatzner fc McKed show that e is quite insensitive to the 
ratio of infall and outflow momentum fluxes. Nevertheless, 
0.5 is only an estimate and e could well vary during accre- 
tion. This is especially true if the protostellar wind were ever 
to truncate accretion, as e{t) — > when this happens. We 
briefly consider other values here. 

The primary effect of varying e, while fixing M*f, E c i, 
and k p , is to change the core mass required to make a star of 
that mass. Suppose we halve e, so that M c must double. The 
accretion time then increases, and M+d decreases, by a factor 
2 1/4 (for k p = 3/2). This mildly stabilizes the disk. But at 
the same time, R c has been increased by 2 1//2 , j has gone 
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Figure 4. Effect of varying the star formation efficiency e in 
the fiducial core collapse model. Dashed lines: critical disk radius 
Rd crit f° r fragmentation; solid lines: expected disk extent R d . 
Intersections are as marked. 




M t (M ) 



Figure 5. Growth of disk and critical radius with the mass of 
a protostar accreting toward 3OA/0 in the fiducial core collapse 
model. In addition to the expected disk radius, we show the 
splashdown radius of the inner infall streamline, calculated as- 
suming e = 0.5. 



up by 2 3 ' 4 , and Rdj has expanded by 2 3//2 . Balancing these 
contributions, we expect lowering e to destabilize the disk. 
This was predicted also in equation (|35[) , where lowering e 
is seen to decrease stability in an active disk. Passive disks 
are extremely insensitive to e. 

Figure @ corroborates our expectation by showing that 
lower values of e correspond to less stable disks. Indeed, 
the mass at which fragmentation sets in is sensitive to e, 
specifically, M cr i t <x e 2 ' 6 , while the critical disk radius is 
relatively constant. Does this mean that a decline in core 
efficiency over time destabilizes disks? Probably not, since 
most of the mass will be accreted at an intermediate value 
of e (see the discussion below equation [18]) . 



4.3 Time Evolution of Disk Fragmentation 

For those disks that do suffer fragmentation, it is useful to 
know whether this happens early or late in accretion and 
how much matter is potentially affected. To address these 
questions we construct the time history of the accretion rate, 
stellar mass, and disk radius for a star with M*/ = 30 Mq 
in the fiducial core collapse model. For this calculation we 
use the lumin o sity h istory of such a star as presented by 
iMcKee fc Tanl (|2003h . The scalings R d (t) oc Af*(i) 1/(3 - fe " ) 
and M* d (t) oc M*(t) (6-2fe ' ,)/(6 ~ 3fc ' ,) permit us to gauge disk 
fragmentation through time. 

Figure (JS| shows the evolution of R d and Rait during 
accretion. In this plot, the relative constancy of the frag- 
mentation radius is due to the enhanced effect of irradi- 
ation at low masses. We find that a star destined to be- 
come 30M© (born of a 60 Mq core) has a disk that crosses 
into the regime of fragmentation when the protostar has ac- 
creted approximately 5.6Mq. The fact that this is slightly 
higher than the critical mass identified earlier is to be ex- 
pected: the inner 5.6M© of a larger core is equivalent to a 
5.6Mq core with a slightly higher column density - specif- 
ically, E c oc (AL/M if )' 1/3 , for k p = 3/2 (cf. equation [33)) . 
The somewhat higher column density implies a smaller and 



somewhat stabler disk, leading to a slightly higher mass scale 
for fragmentation. 

We also show on figure ([5} the radius of the inner- 
most streamline of in falling material from the envelope from 
iTerebev et al.1 <|l984h . a gain assuming an accretion efficiency 
of 50%. This, along with stochastic variations in disk angu- 
lar momentum about its typical value, suggest that accretion 
can coexist with fragmentation so long as the disk is not too 
far beyond the fragmentation threshold; see £|5 .3 1 for more 
discussion. 



4.4 Disk mass and global instability 

At the typical fragmentation radius of 100 — 150 AU, the 
mass scale of a disk with Q = 1 is 



nR d T, d (R d ) 



n ,„ 1/12 / Rd \ 

°- 10e ° 5 ll50AuJ 



1/2 yil/4 



M 



1/4 ■ 

/,30 



(36) 



Globa l instabilities of the disk are triggered by the total disk 
mass (Ad ams et al.|[l989l ; IShu et aLlll990r ), which is larger 
than ■KR 2 d Y, d (R d ) by the factor 2/(2 - fc E ) if E d oc r _fcs 
within R d . There is therefore the possibility t hat the fast an- 
gular momentum transport by these modes (Lauehlin et al. 
1998) suppresses local fragmentation, but we consider it un- 
likely that fragmentation is eliminated by this process. 



5 CONSEQUENCES OF INSTABILITY 
5.1 Fragment Masses 

Once a disk fragments, what objects form? [Goodman fc Tanl 
(2004) determine the initial fragment mass based on the 
wavenumber of the most unstable mode in the disk. The 
corresponding wavelength of this axisymmetric mode is: 



A(r) = 2tt 



c acl 



ttGE' 



(37) 
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where E and c s are both functions of radius within the disk. 
Assuming that the fragment has compara ble dimension az- 
imut hally, the corresponding mass scale is (jGoodman fc Tar] 
|2004) 



10r 



Mb, 



= A E 
Q G 



(38) 



i.e., roughly the amount of mass accreted in 2Q orbits. We 
assume fragmentation only occurs when Q — » 1. We consider 
our rather idealized estimate of Mf rag uncertain by at least a 
factor of two. In reality the initial mass is likely a stochastic 
variable, best determined from numerical simulations (R. 
Rafikov, private communication). 

Once a fragment forms, its growth is controlled by ac- 
cretion of surrounding gas, migration through the disk, and 
collisional or gravitational interaction with other fragments. 
Rather than address these questions in detail, we will draw 
preliminary conclusions by comparing Mf rag to two critical 
scales: the gap opening mass M gap and the isolation mass 
Mi so . A fundamental uncertainty is the state of the gas disk: 
again, we assume Q — 1. 

When the gravitational torques exerted on the disk by 
the fragment exceed viscous torques, a gap opens around 
the fragment. iRafikovl l|2002l ) estimates 



2c 



3fiG 0.043 



0.37^A/ fras 



(39) 



Since this is less than Mf rag , and gets even smaller if the disk 
viscosity goes down and t hus coolin g time goes up relative 
to the critical state in the Gammi c (2001) simulations, we 
expect fragments to open gaps immediately. 

Gap opening slows but does no t necessarily end ac- 
cretion (|Artvmowicz fc Lubowl ["l996T ). A possible limit to 
growth is set by the point at which t he fragment accretes all 
the mas s within its Hill radius (e.g., i Goodman fc Tanll2004l . 
but see lArtvmowicz fc Lubowl ). This defines the isolation 



M iso (r) 



(27r/ w r 2 E) 



3/2 



9M, 



1/2 



(40) 



where fu ~ 3.5. 

Figure ([6| compares the initial, gap opening, and isola- 
tion masses for a range of M*/, evaluated at 7? cr it near the 
end of core accretion (in the fiducial core model). Because 
gap opening slows accretion, we expect the masses of disk- 
born stars to resemble Mf rag more than Mi so . In this case 
they will be low-mass stars of order 0.2-0.4 Mq. 

Fragmentation tends to set in at ~ 100 — 200 AU, as 
we noted in H2.4I However gap-opening fragments should 
be swept inward if disk accretion continues; only a single 
disk mass of accretion is required to bring them in to the 
central object. If they move inward by about a factor of 
30, t hey will be wit hin a few stellar radii at carbon igni- 
tion (|WebbinkNl9"85h . As candidates for mass transfer and 
common-envelope evolution, such objects can serve as reser- 
voirs of matter and angular momentum for the relic of the 
central star's supernova explosion. 
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Figure 6. Our estimate of the initial fragment mass, compared to 
the gap opening mass and isolation mass, at the end of accretion 
in the fiducial core collapse model. We truncate the calculation 
where dust sublimation in the envelope makes the critical radius 
determination uncertain. 



5.2 Observability of Disk-Born Stars 

We have predicted that stars with extended protostellar 
disks will produce low mass (M5-G5) companions. Thus we 
would expect many O stars, and perhaps some early B stars, 
to have multiple coplanar companions at separations of or- 
der < 100-200 AU, depending on the amount of migration 
that occurs following formation. In the closest clusters, for 
example, Orion, this correspondes to an angular separation 
of approximately 0.4" and an apparent bolometric magni- 
tude difference of, at a minimum, ~ 13 magnitudes. This im- 
plies that even in K-band with an AO system such as that of 
VLT, such objects would be difficult to observe (M. Ahmic, 
private communication). Similarly, the combination of AO 
with a coronograph (e.g. the Lyot coronograph on AEOS) 
can provide a dynamic range of up to up to 8 H-band magni- 
tudes at a few hunderd mas, whic h is still too small to detect 
the aforementioned companions (iHinklev et al. I [20061 ). If not 
detectable during the main sequence life of the primary star, 
the presence of such companions might be observable via bi- 
nary interaction once the primary evolves. 



5.3 Disk Efficiency 

If a 3OM0 star will suffer disk fragmentation early in its ac- 
cretion, as we estimated in i]4.3l on the basis of the turbulent 
core model, then we must address how this might impact 
subsequent accretion. One possibility is that gap-opening 
fragments will be swept inward to the central star, in which 
case its final mass will be u naffected. Th is outcome resem- 
bles the scenario outlined bv lLevin] (|2003h for gravitationally 
unsta ble AGN accretion. Alternatively, iTan fc Blackman] 
l|2005l) suggest low-luminosity AGN may be starved of gas 
by fragmentation. 

In the latter scenario, a strongly unstable disk will have 
a low disk efficiency, Ed = M+fM+d- We can estimate Ed in 
the limit that none of the gas entering an unstable region 
of the disk ultimately accretes onto the star. The splash- 
down radius of the innermost streamline was estimated in 
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expression (|12|l . Given that outflow removes matter from the 
inner streamlines, and that fragmentation removes it from 
the outer portions of the disk, the fraction of mass that suc- 
cessfully accretes (after striking the disk) is 



Sd = 1 



1-1- 



Rc 



Rd 



1/2 



(41) 



Of course, this expression is only applicable when it yields 
< Ed < 1, i.e., when the disk is partly but not wholly 
unstable. 

Two complications arise when evaluating equation (|41[) . 
First, the critical radius J? C rit must be calculated using the 
mass accretion rate outside of itself. Second, recall that the 
angular momentum of infalling gas derives from its initial 
turbulent velocity and is likely to vary in direction and mag- 
nitude. The fluctuations of j were estimated by the distri- 
bution fj, which appeared in equation (|25[) . Accounting for 
both of these effects, we find that accretion can continue 
even in an actively unstable disk. For 30Mq, we find that 
the infall streamline remains within the stable disk radius 
through the end of accretion. Even if all the gas entering a 
fragmenting region is consumed, this need not fully starve 
the central object. However, the tendency to fragment be- 
comes much stronger for more massive stars. This is espe- 
cially true for those (M*/ ~ IIOM0, from Q that accrete 
rapidly enough that their disks are destabilized by the drop 
in dust opacity. 

In reality, we expect some gas to accrete through the 
unstable region. A rough upper limit would be to adopt 
the accretion rate for a Toomre-critical (Q = 1) state, and 
assume that any surplus is consumed by fragmentation. In 
this case, the central accretion rate would be limited by the 
temperature of the coldest region of the disk (according to 
eq. |4]) . Numerical simulations will ultimately be required to 
quantify the behaviour of Toomre-unstable disks. 

Furthermore, there are numerous feedback mechanisms 
that have not been addressed. For example, if fragmentation 
halts accretion, this could change the outflow power, the 
shape of the outflow and infall cavity and thus the heating 
of the disk by reprocessed starlight. This interplay will be 
investigated in future work. 



6 DISCUSSION 

Our primary conclusion f i|2.4[l is that massive-protostar 
disks that accrete more slowly than ~ 1.7x 1O -3 M0 yr _1 are 
subject to fragmentation at disk radii beyond about 150 AU. 
This critical radius is set primarily by the viscous heating of 
the disk midplane as it accretes, with reprocessed starlight 
playing an equal or secondary role for most stellar masses. 
As all the disks we consider are optically thick, the critical 
radius depends on the Rosseland opacity law kr(T) within 
dusty disk gas. 

Comp aring to our conserva tive estimate of the disk ra- 
dius in the lMcKee fc Tanl (|2003h model for massive star for- 
mation by the collapse of turbulent cores (32}, we find that 
fragmentation is marginal for stars accreting four to 15 so- 
lar masses; higher-mass stars are increasingly afflicted by 
disk fragmentation. Although the mass at which fragmenta- 
tion sets in is sensitive to our somewhat uncertain fragmen- 
tation criterion ( i|2.1|) and angular momentum calculation 



( i|3.2l and the Appendix), we have been conservative in five 
ways: (1) by adopting a fragmentation te mperature lower 
than that implied by the I Rice et ail (|2005l ) simulations; (2) 
by adopting a low estimate of the specific angular momen- 
tum that determines the disk radius; (3) by adopting a rela- 
tively opaque model for the disk's Rosseland opacity; (4) by 
ignoring the shielding effect of a moderately opaque infall 
envelope, and (5) by adopting a low estimate for the cooling 
time, and a hard equation of state. All of these approxi- 
mations should, if anything, lead us to underestimate the 
prevalence of disk fragmentation. Along with the existence 
of turbulent fluctuations in j (jJU, these points ensure that 
some massive stars above ~ lOM© experience disk fragmen- 
tation. As noted in i|4.4l we cannot rule out the possibility 
that global instabilities flush disk material fast enough to 
suppress fragmentation, but we consider it unlikely that this 
prevents all fragmentation. 

We therefore expect multiple, coplanar, low-mass (M5 
to G5, H5.1|l companions to form around many O (and some 
B) stars. Given initial separations of order 100-200 AU, 
their photospheric emission is not observable with present 
techniques. Disk migration, followed by mass transfer or 
common-envelope evolution, may however make them evi- 
dent as the primary evolves f q5.2p . 

Even when disks fragment, we expect some accretion 
onto the central star - if only because of material that falls 
within the fragmentation radius f £|5.3[) . Although we cannot 
yet quantify the disk efficiency parameter Ed, we expect it 
to be significantly less than unity for those early O stars 
(M*f ^ 50Afo) whose disks are most prone to fragment. 



6.1 Imprint on the initial mass function 

The stabilizing effect of viscous heating is absent in disks 
that accrete more rapidly than 1.7 x lO _3 M0yr _1 , thanks 
to a sharp drop in dust opacity at ~ 1050 K (see Q. As 
this affects stars of mass greater than 87(t aC c/lO 5 yr)M0, or 
about 110 Mq in the fiducial core model, it may be related 
to the cutoff of the initial mass function (IMF) at about 120 

Mq. 

Several other explanations for this cutoff have been pro- 
posed, all involving the increasing bolometric luminosity, 
ionizing luminosity, or outflow force emitted by the central 
star. Starvation by disk fragmentation has the distinctive 
feature that it becomes much more severe at a specific accre- 
tion rate. For this reason we expect it to produce a sharper 
IMF cutoff. (The transition to super-E ddington luminosities 
could also produce a sharp cutoff, but lKrumholz et al.ll2005l 
argue that this can be overcome by asymmetric radiation 
transfer.) 

Note also that disk accretion is destabilized by rapid 
accretion, whereas rapid accretion quenches the effects 
of direct photon force and of the ionizing radiation 
jWolfire fc Cassinelli|[l987h . Disk fragmentation may close 
an avenue by which very massive stars would otherwise form. 



6.2 Proto-binary disks 

IPinsonneault fc Stanekl (2006) state that close massive bi- 
naries < 10 year periods) are more likely to have nearly 
equal masses. More generally, the binarity fraction overall 



© 0000 RAS, MNRAS 000, 000-000 



12 Kratter & Matzner 



among massive stars is higher than thei r low mass counter- 
parts (1.5 versus 0.5, iBallv et al] <|2005l >). Due to the high 
fraction of roughly equal mass binaires, their effect on disk 
dynamics must be addressed in future work. Moreover, the 
stellar densities in regions of HMSF are high, suggesting that 
other cluster stars might be close e nough to interfer e with 
disks stretching out to ~ 100 AU jBate et alj|2003f ). It is 
not currently clear whether disk accretion can preferentially 
grow a low-mas s companion until its mass riv als that of the 
primary star, as lArtvmowicz fc Lubowl (|l996l ) suggest. If so, 
then disk fragmentation may be relevant in the production 
of equal-mass binaries; if not, they must form by another 
mechanism. In any case, the multiplicity of the center of 
gravity should be accounted for in future work. It seems 
unlikely, though, that a binary with ~ 10 year period will 
stabilize the fragmenting regions whose periods are ^ 300 
years. 
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Profile 


/8 = l/4 


(3=1/2 


turbulent core 


0.2704 


0.4206 


critical Bonnor-Ebert sphere 


0.2730 


0.3828 


uniform region (k p = 0) 


0.3430 


0.4714 


thin shell 


0.4714 


0.6324 



Table 3. Values of ij 2 ) 1 ^ 2 /(R(c 2 } 1 ' 2 ) in our model for turbulent angular momentum. 
7 APPENDIX 

We present here an estimate of the angular momentum of cores initially supported by turbulent motions, generalizing the 
results of ML05 to arbitrary line width-size relations. The analytical treatment of this problem rests on several idealizations 
about turbulent velocities: (1) that they are isotropic and homogeneous; and (2) that the Cartesian velocity components are 
neither correlated with each other, nor (to any appreciable degree) with density fluctuations. For simplicity we also assume 
(3) that the core density profile is spherically symmetric and can be captured in a single function p(r). When evaluating our 
formulae we assume (4) that the velocity difference between two points scales as a power of their separation, and (5), for the 
purpose of computing fluctuations, that the velocity components are Gaussian random fields. 

One might object that condition (1) is inconsistent with the turbulent support of an inhomogeneous density profile. 
Consider, however, that if the core profile is a power law p(r) cx r~ kp , then the turbulent line width must scale as cr{r) tx r 13 
where 



/3 = 1- k p 12. 

A velocity field with this scaling is consistent with our conditions (1), (2), and (4) if 
([v I (r 1 )-v,(r 2 )] 2 ) = fc|r 1 



- r 2 \ 2l3 Sij 



(42) 



(43) 



where angle brackets represent an ensemble average and k is a normalization constant related to the virial parameter a. ML05 
considered the case (5=1/2 appropriate for giant molecular clouds (k p — 1); we generalize their formulae to other values of 
f3, including the fiducial value j3 = 1/4 corresponding to k p = 3/2. 

Our goal is to compute the expectation value and dispersion of the core specific angular momentum j, normalized to R C <J 
where a is the one-dimensional line width of the core. We find that, under our assumptions, 



(f) = - 



d 3 np(n) d 3 r 2 p(r 2 ) 



M 



Al 



ri ■ r 2 ([u a; (ri) - v x (v 2 )\ ) 



(44) 



and that 

(cr ) = — 

x 1 2 



d 3 np(n) d 3 r 2 p(r 2 ) 



M 



M 



(Mn) -v x {r 2 )f) 



(45) 



where M — J pd 3 r and all integrals are restricted to the region of interest (typically the core interior) . These equations, which 
we will prove below, agree with the formulae in ML05's Appendix but allow ({v x (ti) — v x (r 2 )] 2 } to be an arbitrary function of 
ri — r 2 |. ft is important to note that the two formulae are identical up to the factor — 2ri • r 2 in the integrand. The negative 
sign in equation (|44p ensures that (j 2 ) is positive, since ([^(ri) — v x (r 2 )] 2 ) takes higher values when |ri — r 2 | is large, hence 
when ri • r 2 is negative. 

To evaluate equations (|44p and ((45} we make use of spherical symmetry and impose equation (|43[1 for the velocity 
correlations. Defining fj, = ri • r 2 /(rir 2 ) as the cosine of the angle between ri and r 2 , 



if) 

and 

<- 2 > 



k / dri47rr 1 



dri 



p(n) 
M 



j 2 p{r 2 ) 
dr 2 2irr 2 — 



r\r 2 ir 



r+r-i) "J" 



dp(l - qpf p 



2-K 2 k 
M 2 



dr 2 (rl +r 2 ) 2+l3 nr 2 p(n)p(r 2 ) 



(1 + 9) 



(3+1 



(1-9) 



(3+2 



(1 + 9) 



0+1 



(1-9) 



13+1 



2 

2n 2 k 
M 2 



M 



j r, 2p(r 2 ) 
dr 2 2nr 2 — 



2\/3 
r 2 



(3 + 2 



dp(l -qp,y 



(3 + 1 



(46) 



dri 



dr 2 (r\ +- r 2 



i+0 



rir 2 p(ri)p(r 2 ) 



(l + qf +1 -(l~ q r + i 



(3 + 1 



(47) 



where q = 2r\r 2 j(r\ + r 2 ), as in ML05. Table[7]and figure [7] provide values of (j 2 ) 1 ^ 2 / ^(i?(<r 2 ) 1//2 ) for several density profiles, 
both for the internal spectrum given by equation (|42|) . and for (3 = 1/2, which may better represent the background spectrum. 

t igure □ plots the ratio (j 2 ) 1/2 /(R(a 2 ) 1/2 ) as given by equations (|46|l and (|47l) for a turbulent core profile, a critical 
Bonnor-Ebert sphere, and a thin shell. The results are very close to power laws in (3: 
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0.1 0.2 0.3 0.4 0.5 

Line width-size exponent 



Figure 7. Values of (j 2 ) 1 / 2 /(/(.(cr 2 ) 1 / 2 ) evaluated for turbulence with line width-size exponent fj and three relevant density profiles. For 
hydrostatic turbulent cores, p oc r~ 2 ( 1— W, 



2\l/2 



0.655/3 0,638 singular turbulent core 
0.537/3 ' 488 critical Bonnor - Ebert sphere 



J?{a2>i/2 " I 0.648/3 u ' 4bH uniform - density sphere 



(48) 



0.849/3 0424 thin shell. 



The spin parameter 6j makes reference to the velocity dispersion a(R) at radius R, rather than the mean velocity dispersion 
(cr 2 ) 1/2 within R. For this we use the scaling a 2 oc GM/r oc M 2l3/(3 ~ kp) to derive 

a (R) 2 -i , 2/3 

/ OX ^ "T" 



3 k 



With this correction factor accounted for, 9j / fj is still very close to a power law of j3: 
( 0.504/3 0,552 singular turbulent core 



j- - { 0.987/3° 651 uniform - density region 
h { 0.849/3 424 thin shell. 



(49) 



(The last line is unchanged, as there is no correction for a thin shell.) These formulae were quoted in £13.21 

We now return to the derivation of equation (|44|l ; equation (|45[) is treated below. The z component of the specific angular 
momentum of the core is 



J 2 



d d ri / d 3 V2^{x 1 v yl - yiv xl )^-(x 2 v y2 - y 2 v x2 ) 



Pi , 



' M 



M 



^2 / d 3 ri / d 3 r 2 pip 2 (x 1 x 2 v y iv y2 + yiy 2 v x iv x2 - xiy 2 v y iv x2 ~ yix 2 v x iv y2 ) 



(50) 



here, as below, we use subscripts 1 and 2 to indicate functions of ri and r 2 , i.e., pi = p(ri). On ensemble averaging of the 
second line, the last two terms in equation (|50p become zero thanks to the <5y in equation (|43[1 ■ The first two terms are equal, 
thanks to spherical symmetry; thus 



(jz) = -pj / d3ri / d3r 2PiP23:ix 2 (v yl Vy 2 ). 



(51) 



Spherical symmetry allows us to replace xix 2 with n • r 2 /3 and j 2 with j 2 /3, so 
0' 2 > = J~p J d3ri J d 3 v 2 p 1 p 2 Y 1 -r 2 (v yl Vy 2 ) 

— j d 3 ri / d 3 r 2 pip 2 ri ■ r 2 ((u 2 !) + (t; 2 2 ) - {(v yl - w !/2 ) 2 )) . 



(52) 



The first two terms within the brackets are zero, as spherical symmetry requires them to be even in the space coordinates 
while ri • r 2 is odd; the surviving term yields equation (I44p . 
To prove equation (I45|l we write 

Ma 2 = j d 3 rp{r)(v x - v x f = / d 3 rp(r)(« 2 - v 2 x ) (53) 
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where v = J vp(r)/Md 3 r is the center of mass velocity. In equation (|50[1 we replace v% with J j vi ■V2(pi/M)(p2/M)d 3 Tid 3 r2, 
and we bring the first term into similar form by multiplying it by J (p 2 / 'M)d 3 Y2 = 1: 

r 

d 3 ri / d 3 r 2 pip2^i(^i - ^2) 



d ri / d r 2 piP2(v x i - v x2 ) . (54) 



The second line follows from the first by noting that the integral is antisymmetric under interchange of ri and r 2 . When 
averaged, it yields equation (J45J . 
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